Time-dependent electron transport through a strongly correlated quantum dot: 
multiple-probe open boundary conditions approach 
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We present a time-dependent study of electron transport through a strongly correlated quan- 
tum dot. The time-dependent current is obtained with the multiple-probe battery method, while 
adiabatic lattice density functional theory in the Bethe ansatz local-density approximation to the 
Hubbard model describes the dot electronic structure. We show that for a certain range of voltages 
the quantum dot can be driven into a dynamical state characterized by regular current oscillations. 
This is a manifestation of a recently proposed dynamical picture of Coulomb blockade. Further- 
more, we investigate how the various approximations to the electron-electron interaction affect the 
line-shapes of the Coulomb peaks and the I-V characteristics. We show that the presence of the 
derivative discontinuity in the approximate exchange-correlation potential leads to significantly dif- 
ferent results compared to those obtained at the simpler Hartree level of description. In particular, 
a negative differential conductance (NDC) in the I-V characteristics is observed at large bias volt- 
ages and large Coulomb interaction strengths. We demonstrate that such NDC originates from the 
combined effect of electron-electron interaction in the dot and the finite bandwidth of the electrodes. 

PACS numbers: 05.60.Gg, 71.10.Fd, 73.23.Hk 
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I. INTRODUCTION 

Electron transport through nanoscale devices is a di- 
verse subject, which is currently the focus of extensive 
experimental and theoretical research. The fuel of such 
interest is the expectation that nanoscale objects, such 
as quantum dotsi and even single molecules j£ are to 
become active components in novel electronic devices, 
which potentially offer unique advantages over existing 
technologies j 3 - At the fundamental level, the physics of 
such reduced-dimensional systems is dominated by quan- 
tum effects. Among them are electron correlations, which 
strongly affect the electron transport at this level of con- 
finement, giving rise to prototypical quantum phenom- 
ena, such as Coulomb blockade^ 5 - and the Kondo ef- 
fect^ 

While the Landauer formula is the solution to the non- 
interacting quantum transport problem^ the interact- 
ing case continues to be challenging to the theory. The 
latter is typically approached with the non-equilibrium 
Green's function (NEGF) formalism^ which allows, in 
principle, the derivation of an interacting many-body 
Landauer-type formula for the steady-state current in 
the case where interaction is limited to a finite region 
in spaced In practice, for the majority of the state- 
of-the-art ab initio transport calculations and numeri- 
cal algorithms ) 12 " 15 the method of choice for the elec- 
tronic structure description is the density functional the- 
ory (DFT). However, typical steady-state DFT+NEGF 
transport schemes have a range of limitations, both con- 
ceptual and technical^ 

At the fundamental level it has been recently demon- 
strated, at least for the case of a single Anderson impu- 
rity model, that the linear response conductance calcu- 
lated from the Kohn-Sham levels for the exact exchange- 
correlation (XC) functional reproduces closely that com- 



puted with many-body approaches! 17 : 18 If the same holds 
true for ab initio DFT, then the DFT+NEGF scheme will 
provide a complete solution for the zero-bias limit. Still, 
on the practical side, the commonly used approximations 
to the XC functional, lacking the so-important derivative 
discontinuity^ fail to capture essential physics for the 
transport in molecular junctions, qualitatively mispre- 
dicting the conduction rcgimei 20 : 21 Different is the situ- 
ation at finite bias, where, let alone the implementation, 
conceptual concerns reflect on the very applicability of 
a ground-state electronic structure theory to an intrinsi- 
cally non-equilibrium problem especially if electron cor- 
relations are significant ! 22 ' 23 

One strategy to avoid some of the shortcomings of 
using equilibrium DFT has been sought in its natural 
extension, time-dependent (TD) DFT^ with practical 
schemes for TD transport having been developed^ In 
general, real-time TD schemes for quantum transport can 
be roughly divided into two categories based on their as- 
sumption for the initial conditions. In one case the elec- 
trodes arc prepared in equilibrium with the poles of a 
battery, but not yet connected to the nanoscopic device. 
The current then starts to flow when the connection is 
made. In the other the system electrodes+device is ini- 
tially at equilibrium and subsequently an electric field is 
applied to the electrodes. The former assumption, where 
two initial electrochemical potentials are well defined, is 
more in the spirit of the Landauer transport picture. The 
latter is instead more DFT-friendly, as the starting point 
is the ground state of the system^ 

There has been evidence that these two TD trans- 
port variants agree in the non-interacting case, i.e. they 
lead to the same history-independent steady-state cur- 
rent More recently, the latter variant combined with 
the TDDFT, further equipped with a novel XC functional 
carrying the physical derivative discontinuity, has been 
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applied to study the transport through a quantum dot 
in the Coulomb blockade (CB) regime by Kurth et al. 
in Ref. [I?]]. In particular that work has put forward an 
important novel description of CB as a dynamical pro- 
cess with rapidly oscillating local currents, inaccessible 
by conventional steady-state transport models. 

In this work we adopt another recently proposed 
TD transport scheme, the so-called, multiple-probe bat- 
tery (MPB) method ; 28 ' 29 to study electron transport 
through a strongly correlated quantum dot. The MPB 
scheme was first proposed in the context of corre- 
lated electron-ion dynamics and was applied to a wide 
range of problems, such as current-induced heating in 
atomic wiresi 28 ' 30 This method belongs to the first of 
the fore-mentioned categories and enables the realization 
of an external battery within the finite system of elec- 
trodes+device. The external bias is introduced through 
the difference in the electrochemical potentials of the set 
of reservoirs, or probes, attached individually to each 
atom in a pair of large but finite metallic electrodes 
(leads). The scheme is very tractable computationally 
and has the control knobs to be an arbitrarily close ap- 
proximation to the non-interacting Landauer transport 
in the long-time dc limit. 

The MPB time-propagation scheme is based on the in- 
tegration of the Liouville-von Neumann equation of mo- 
tion for the reduced density matrix of the system, in 
which the open boundaries are described explicitly by a 
source and a drain term. For the TD Hamiltonian of the 
quantum dot, entering the equation of motion, we adopt 
the description used by Kurth et al.^2. This is based on 
the adiabatic Bethe ansatz local-density approximation 3 ^ 
(adiabatic BALDA, or ABALDA) to the XC functional, 
which exhibits a derivative discontinuity at half-filling. 

By investigating the real-time evolution of the current 
through the quantum dot, we find an agreement with 
Ref. |27| , i.e. for a certain set of parameters the system 
does not reach a steady state but rather remains in a 
dynamical state, characterized by oscillations in the cur- 
rent. Furthermore, we try to interpret the TD results in 
terms of the more familiar steady-state picture of trans- 
port. In particular, we construct the current- voltage, I- 
V, characteristics of the quantum dot from the long-time 
average of the current and the voltage obtained from the 
TD simulations. This is done for a wide range of pa- 
rameters, even in the cases when a steady state is not 
achieved. Importantly, we observe a drop of the current 
as a function of the source-drain voltage and, as a conse- 
quence, a negative differential conductance (NDC) above 
a critical bias voltage. We demonstrate that such an ef- 
fect is not possible if the derivative-discontinuity is not 
included in the one-particle potential. 

This is particularly interesting in view of some recent 
contrasting results. On the one hand a number of stud- 
ies, based on several distinct many-body approaches^ - — 
attribute the NDC mainly to electron-electron interac- 
tion. On the other hand, it has been demonstrated by 
Baldea and Koppel 3 -^ that in the case of an exactly solv- 



able model for a non-interacting dot within the steady- 
state formalism, the finite bandwidth of the electrodes 
can alone lead to pronounced NDC for a wide range of 
parameters. Here we find a numerical proof that this re- 
sult can be generalized to the interacting case and time- 
dependent transport. Our calculations suggest, however, 
that for the system considered here, the NDC is due to 
a combination of two effects, namely electron-electron 
interaction on the dot and the finite bandwidth of the 
electrodes. 

Our paper is organized as follows. In the next sec- 
tion we introduce the model system and our theoretical 
framework, i.e. the Hamiltonian and the computational 
scheme for MPB quantum transport. In the first part of 
Section IIIII the I- V characteristics of a non-interacting 
quantum dot calculated by using the TD-MPB method 
is compared to analytic NEGF results. We then discuss 
the finite electrode bandwidth as a source of NDC. In the 
second part of Section IIIII we present the TD results for 
a strongly correlated dot in the CB regime. Finally, we 
propose an explanation for the observed NDC in the I- V 
characteristics. 

II. METHODS 

The model system considered in this work is presented 
in Fig. [TJ This consists of a central region, which con- 
tains the quantum dot surrounded by two Ad-site long 
atomic chains at both sides, and two one-dimensional fi- 
nite leads, each counting N-^m-j atoms. The physics of 
the quantum dot connected to two leads is described by 
the Anderson impurity modeh 11 ' 36 The Hamiltonian of 
the total system thus reads 

Hs= H a + H T + H QD . (I) 

a— L.R 

Here the first term is the nearest-neighbors single-orbital 
tight-binding (TB) Hamiltonian describing respectively 
the left-hand side (a=L) and right-hand side (a=R) lead. 
This is written as 

H a =Yl Sia d Cf a + 70 ( Cl £ i+l a + h - c ) , ( 2 ) 
i,(T i,a 

where Ei a are the on-site energies and 70 is the hopping 
integral; c^j(cf Q ) is the creation (annihilation) operator 
for an electron with spin a (<7=t, i) at the atomic site 
i of the lead a (the index i = 1, ..,N a runs from left to 
right for q=R and from right to left for a=L). Note that 
two atomic chains on each side of the quantum dot are 
also described by a TB model with the hopping integral 
70 and therefore they are included in the Hamiltonian of 
the leads. 

The second term in Eq. ((T|) describes the tunneling 
between the quantum dot and the two adjacent sites and 
it is given by 

Ht = Y^( $ £ il + $ %r + h.c) , (3) 



3 



QD 



sd 



S S + + 4 4 



I 



FIG. 1: (Color online) Schematic of the model system consid- 
ered in this work: the central region consists of a quantum dot 
(QD) surrounded by two iVd-site long atomic chains, which in 
turns are attached to two one-dimensional leads comprising 
respectively JVl and JVr sites. Here 70 is the hopping integral 
in the leads and in the two chains, and 7 C is the lead to dot 
hopping. V g denotes the gate voltage, acting locally on the 
dot, and V s( j is the source-drain voltage applied across the 
entire system. 



where Cq (cq) is the creation (annihilation) operator for 
an electron with spin a on the dot and ~f c is the hopping 
integral between the dot and site i=l in the lead a. 
Finally, the Hamiltonian of the quantum dot reads 



(4) 



where V g is the on-site energy of the dot, which acts as 
a local gate voltage; U (U > 0) is the charging energy, 
which expresses the strength of the Coulomb repulsion 
on the dot; Uq =c^ Cq is the site-occupation operator. 

Within the lattice DFT framework^ the many-body 
Hamiltonian in Eq. (U) is mapped onto an effective single- 
particle Kohn-Sham Hamiltonian which, in the local den- 
sity approximation, reads 



vks [nol n . 



(5) 



Here no is the charge density of the dot and vks is the 
effective Kohn-Sham potential, which can be written as 
a sum of three terms 



«ks [no] = V e + -£-U + vxc [no 



(6) 



The second and third terms are respectively the Hartree 
and the XC potential. The latter is approximated by 
a modified BALDA potential, specifically tailored to 
a nonuniform configuration with a weakly coupled dot 
(we refer to Ref. [27| for the exact expression and the 
parametrization) . 

Notably, such vxc exhibits a derivative discontinuity 
at no=l, i.e. at the phase transition of the ID Hubbard 
model. In practice, however, we use a continuous approx- 
imation to the BALDA potential 2 ^ where the true dis- 
continuity, expressed through a Hcaviside step function 
6(no), is replaced by a function /(no) = l/(e^ rao_1 ^ a + 1) 
with a being a smoothing parameter. We use a = 10~ 7 , 
which guarantees a very sharp slope at no = 1. In 
our simulations we consider three levels of description: 
(i) U = 0, or non- interacting case, for which the effec- 
tive potential of the dot is simply given by VKS=V g , (ii) 



vxc ^ 0, or the Hartree approximation, where the po- 
tential on the dot is vn=V g + Uuq/2; and (Hi) the full 
discontinuous effective potential, given by Eq. (|6]), which 
we refer to as i>ks for clarity. 

In order to introduce the time-dependence in the 
Hamiltonian of the quantum dot, we use the adiabatic 
approximation, where vq is assumed to depend on time 
only through the instantaneous charge density of the dot 



VKs(t) = v K s[no(t)] 



(7) 



The question of the applicability of such adiabatic lo- 
cal approximation to the description of non-equilibrium 
transport in strongly correlated systems has been ad- 
dressed in recent two works respectively by Uimonen 
et al£^ and Khorsavi et al£L In particular, a compara- 
tive study between the TDDFT approach with ABALDA 
(TDDFT+ABALDA) and the many-body perturbation 
theory, applied to out-of-equilibrium Anderson impurity 
model, has been carried out in Ref. [H)]. The results ob- 
tained with both approaches have been tested against nu- 
merically exact results produced by time-dependent den- 
sity matrix renormalization group theory. It was found 
that, in general, the TDDFT+ABALDA approach is in 
good qualitative agreement with many-body perturba- 
tion theory over a wide range of parameters. However, 
in many cases it overestimates the steady-state currents. 
This problem was linked to the shortcomings of the local 
approximation to the XC functional and, in particular, to 
the absence of electron correlations inside the electrodes. 
Moreover, it was demonstrated in Ref. [39[ that the inclu- 
sion of dynamical correlations, or memory effects, might 
eliminate the multistability in the density and the cur- 
rent, which can be found within the TDDFT+ABALDA 
approach. These are strong indications that more ad- 
vanced non-local, both in space and time, approxima- 
tions to the XC functional are required. However, as was 
demonstrated in Ref. [I?} and as it will be shown in this 
paper, the ABALDA already provides valuable insights 
into time-dependent transport in strongly correlated sys- 
tems. 

We now discuss, following the work of Todorov and 
co- workers ; 28 ' 29 how the open boundary conditions are 
introduced in the MPB setup. In the MPB method, each 
atom i of the leads (with the exception of the atoms at 
both sides of the quantum dot) is connected to an exter- 
nal probe Pi (see Fig. [IJ . All the probes attached to the 
sites in the left (right) lead are kept at the electrochem- 
ical potential /il (mr) an< 4 are occupied according to the 
Fermi-Dirac distribution /l (/r). The source-drain volt- 
age V s d is introduced as V s d = u^ — (here V sr \ is in units 
of eV). For symmetrically applied bias /^l = £f + K*d/2 
and nil = £f — Kid/2, where £f is the Fermi level of 
the electrodes (assumed identical). The time-dependent 
equation of motion for the density matrix of the system 
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coupled to the probes reads 



ihps(t) 



H s (t),ps(t) +t+ p s (t)- p s (t)XT + (8) 



dE . 



The last two terms on the right-hand side are extraction 
(drain) and injection (source) terms, respectively; G + 
(G~) is the retarded (advanced) Green's function of the 
system and it is given by 



G 



Eh-H, n -E ± ±iI,A 



(9) 



where H Sa = J2 a =L,R H a + H T + J2 a V s n% is the time- 
independent part of Hg(t) and A is a dephasing factor 
(see later for an exact definition). The self-energies due to 
the presence of the external probes and the in-scattering 
self-energy are written as 



. r f . r . 



(10) 



^h(E)i L + ^f R (E)i R , (ii) 

Z7T Z7T 



with the broadening T defined as T = 2ir r ypd, where 7p 
is the coupling to the probes, assumed to be identical 
for all sites in the leads, and d is an energy-independent 
constant, which represents the surface density of states of 
the probes within the wide-band limit; Im is the identity 
operator in region M (M=L, R, S). 

Equation (J8J) is derived from a general Liouville-von 
Neumann equation for the total density matrix of the sys- 
tem and the probes combined. It incorporates two main 
approximations: (i) the wide-band limit in the probes 
and (ii) the decoherence in the injection process, intro- 
duced through the relaxation time ta, with A = h/r^ 
[see Eq. ©]. The second approximation essentially de- 
couples, over the time interval ta, the injection of elec- 
trons from the probes into the leads and their subsequent 
scattering from the time-dependent potential inside the 
central region, provided that the latter is long enough. In 
other words the dephasing factor imposes a restriction on 
the size of the central region (2Nd + 1 sites). Therefore 
the inclusion of buffer sites on both sides of the dot 
is essential within the time-dependent formalism. 

The value of A is determined in such way that the dis- 
tance traveled by the electrons during the time interval 
ta is smaller than the distance between the electrodes 
and the interior of the central region, i.e. the quantum 
dot. This condition can be written as v e TA < N^a, where 
V e is the electron group velocity and a the lattice con- 
stant (a = 1). In practical terms, the introduction of 
the dephasing factor allows one to write down the in- 
jection term, which is in general non-local in time, in a 
rather simple time-independent form [see Eq. ([5])]. This, 
however, also introduces an additional broadening, pro- 
portional to A, in the steady-state I-V characteristics, 



which is absent in the standard static NEGF formalism. 
We note that in the steady-state MPB formalism, the 
NEGF result is recovered in the limit of infinitely long 
leads and weak lead-probe coupling^ 

In order to investigate the open-boundary electron 
dynamics in the time domain, Eq. (O is numerically- 
integrated using the fourth-order Runge-Kutta (RK4) al- 
gorithm!^ As initial condition, we use the density matrix 
Ps(to) of an isolated system (not coupled to the probes), 
constructed from the cigenstates of the Hamiltonian Hs ■ 
The open boundary terms are switched on over a short 
time interval of 5 fs and maintained throughout the simu- 
lation. The current through the dot is then calculated as 
a bond current between the dot and the adjacent site4i 
The typical parameters of the MPB setup used in our 
simulations, unless specified otherwise, are -/V L / R = 90 
and iVd = 20. We have tested that further increasing the 
size of the system does not lead to significant difference in 
the I- V characteristics. In order to have one free param- 
eter instead of two, we use the condition A=r/2, which 
has been discussed in detail in Ref. [28|, and T = 0.35 eV 
in our simulations. 



III. RESULTS 

A. Non-interacting case 

As a test of the applicability of the TD MPB method 
we first examine the non-interacting case (U = 0). For 
this situation, we directly compare the I-V characteris- 
tics obtained from the time-dependent simulations to the 
ones calculated by using the standard NEGF-bascd Lan- 
dauer solution, which we refer to as exact NEGF^ The 
comparison is presented in Fig. [5J where the current is 
plotted as a function of the source-drain voltage for the 
non-interacting level aligned with the Fermi level in the 
leads (V g = 0). In the case of the TD MPB approach, 
the value for the steady-state current is obtained from the 
time-dependent simulation for the corresponding value of 
V s d after the steady-state has been established, i.e. when 
the variation of the current with time becomes negligi- 
ble. In the case of the exact NEGF method, we use the 
well-known analytical expression for the non-equilibrium 
current through a non-interacting resonant level coupled 
to two semi-infinite electrode a 10 ' 35 



-Ten = 



2e 



dE 



T L (E)T R (E) 



[E-V s 



-ME)Y 
ME)}. 



[T(E)/2] 



(12) 



Here A(E) = A L {E)+A R (E) and T{E) = T L (E)+T R {E) 
represent, respectively, the real and imaginary part of the 
total self-energy due to the presence of electrodes, with 
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A L (R) and r L(R) given by 



A L (R)(£) = ^|£l(R), 



(13) 



r L(R)(£) = ^^(270 - |^L(R)|)^4 7 o 2 - El {R) ,(14) 

where E a = E — e a , e a being the on-site energy in the 
lead (a =L, R). 



v = o 




exact NEGF 7- V. The agreement is particularly good in 
the highly conducting limit. The smearing of the abrupt 
I-V features at low bias and again the NDC drop at 
V s d < 47o for the weakly conducting limit are inherent 
to the TD MPB method^ These are due to the explicit 
dephasing factor, which simplifies the equation of mo- 
tion for the density matrix by eliminating temporal non- 
localities of the injection. 

In order to eliminate the drop in the current at large 
bias voltages and to focus on the electron interaction at 
the quantum dot, we will use the £l(r)=0 limit in all the 
further calculations presented. In this case, the satura- 
tion current at high voltages is entirely determined by 
the position of the resonant level, set by the gate voltage 
V E (see the inset of Fig (2), relatively to the electrodes 
band center. As the resonant the level approaches the 
band-edge of the leads (V g < 270), the saturation cur- 
rent decreases. In Section IIII B 21 we will recognize the 
contribution of the latter effect to the drop in the cur- 
rent as a function of the source-drain voltage. 



B. Interacting case 



Kh (eV) 



1. Time- dependent transport 



FIG. 2: (Color online) Current through the quantum dot, 7o, 
as a function of the source-drain voltage, V B d, for zero gate 
voltage (V g =0), calculated using the both exact NEGF and 
the TD MPB method, and for two configurations of the leads: 
£ L/ /r,=0 an d £L/R=±F s d/2. The inset shows the I-V charac- 
teristics obtained with the TD MPB approach for e L / R =0 
and three different gate voltages, F g =1.0, 1.5 and 2.0 eV. 
The following parameters are used: 70=— 1.0 eV, 7c=— 0.1 eV 
and £f=0. The source-drain voltage is applied symmetrically, 
/i L/ / R =£F±14d/2. In order to achieve a better agreement with 
the exact NEGF results we use the improved MPB setup with 
iV L(R) = 250, iV d = 90 and T = 0.15. 

We consider two possible limits for the on-site energies 
in the electrodes: (i) the highly conducting regime with 
e Q =0 for all atoms in a =L, R and (ii) the weakly con- 
ducting regime for which the on-site energies in L(R) arc 
shifted in accordance with the respective electrochemi- 
cal potential, £L(R)=±K;d/2. As expected, the difference 
between the I-V curves calculated in these two limits 
becomes significant at large bias, since the transmission 
in case (ii) rapidly drops to zero once the bias voltage 
exceeds the bandwidth of the leads (4|7o|). This high- 
bias NDC effect, stemming entirely from the finite elec- 
trode band-with, is a well-understood feature of steady- 
state transport in low-dimensional yet uncorrclated elec- 
tron systems.— We also note that the low-bias agree- 
ment between the two transport limits can, in principle, 
be extended to arbitrarily high biases V s d by increasing 
70 > V sd /4& 

An encouraging result is that for both the transport 
limits the TD MPB method reconstructs rather well the 



While in the non-interacting case the TD current 
through the dot always reaches the steady-state, in the 
case when electron-electron interaction is considered this 
is not guaranteed. In fact for certain values of the source- 
drain voltage, for which the charge density of the dot 
approaches unity, the system is driven into a dynamical 
state, where current, density and on-site potential oscil- 
late^ without ever reaching a steady-state. 

The question we address here is whether such dynam- 
ical state can be captured by the MPB method. The 
results of our calculations are shown in Fig. [3J For all 
values of the source-drain voltage below a critical value 
V£[ a steady-state is achieved. However, for source-drain 
voltages above V£[, oscillations indeed develop in all 
transport-related quantities. As shown in Fig. [3] for this 
range of V s d the density quickly reaches a critical value of 
no = 1. At the same time the first jump of the on-site po- 
tential occurs, followed by a series of almost rectangular 
pulses [see Fig. [21b)] . Due to the derivative discontinu- 
ity at no = 1 , the on-site potential reaches an oscillating 
regime, abruptly alternating in time between two values, 
one just below and the other just above the discontinu- 
ity. This translates into oscillations of the charge density 
around no = 1 [see the inset in see Fig. 02a)] and also 
into oscillations in the current [Fig. EJb)]. 

Below, we elaborate on the dynamical features ob- 
served for different values of V s d ■ The height of the pulses 
in VKs(t) is equal to the height of the jump of wks[«o] 
at the derivative discontinuity and it is mainly governed 
by the value of the charging energy U. The width of 
the pulses increases with increasing V s d- This essentially 
means that for larger V s d the system tends to stay longer 
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FIG. 3: (Color online) Real-time evolution of the quantum 
dot: (a) Charge density of the dot (no) for four different values 
of the source-drain voltage, V s d=1.3, 1.6, 1.7, and 1.9 eV. The 
inset shows the fluctuation of the density around unity, Sno, 
defined as Sno = (no — 1) X 10 3 . (b) Current through the dot, 
Jo, for two values of V a d'- Vsd=l-6 eV (black solid line), which 
corresponds to the oscillating regime, and V s d = l-3 eV (black 
dashed line) where no oscillations are observed. Note that 
the corresponding Kohn-Sham potential (vks) [red solid line] 
is also in the oscillating regime (V£d=1.6 eV). The following 
parameters are used: 70= — 1.5 eV, 7c=— 0.3 eV, £f=1.5 eV, 
U=2.0 eV, £l(r)=0 is taken as a reference of energy. The 
source-drain voltage is applied asymmetrically (/^l=£f + V s d, 
ur=£f). 



in the state with a larger on-site potential, corresponding 
to the density above 1 . Further increasing V^d will finally 
lead to a steady-state. The exact value of the threshold 
voltage, V£[, is difficult to determine since the on-site po- 
tential changes with time. From simple considerations, 
however, we established that VC" > 1>ks[?\], where n is a 
value of the charge density just below 1. For the set of 
parameters used here V^f « 1.5 eV. 

As discussed by Kurth et al., the dynamical state of the 
quantum dot described above is a manifestation of dy- 
namical Coulomb blockade. By applying a large enough 
source-drain voltage the dot can be charged. However, 
when the charge reaches the critical value uq = 1 , the on- 
site potential immediately increases by an amount, deter- 
mined by Coulomb repulsion U, thus preventing further 
charging. This essentially corresponds to the CB regime. 
In addition, the time-dependent simulations reveal that 
in this regime the quantum dot is alternating between 
two states, separated by an energy barrier determined 
by U. These two states correspond to the fluctuation of 
the charge on the dot around no = 1, which originates 
from the fact that the ABALDA potential has a deriva- 
tive discontinuity at no = 1 but it is a smoothly varying 
function of uq away from this occupation. 

It follows from the discussion that the dynamics of the 
quantum dot in the CB regime, calculated with the TD 
MPB method, is in a good agreement with the results re- 
ported in Ref. p7| both qualitatively and quantitatively. 



We have established numerically that the two different 
methods reproduce practically identical dynamical tra- 
jectories for all the observables in the long-time limit in 
the case of an interacting system. The remaining differ- 
ences arc limited to the early stage of the time-evolution. 
A characteristic feature of the on-site potential of the dot, 
observed in Ref. [27J , is a transition period just after the 
start of the oscillations, where the series of rectangular 
pulses in the time-dependent vks is preceded by a larger 
pulse whose width increases with V^d- This characteris- 
tic transient pulse is not present in our calculations (see 
Fig. E). 

In order to establish to what extent the transient pulse 
is determined by the initial conditions, we performed TD 
simulations for the same system as shown in Fig. [1] but 
without attaching the external probes, i.e. for a closed- 
boundary finite system. Instead, we applied the source- 
drain voltage as a rigid shift of the on-site energies in the 
left lead, i.e. a term V s d J2i a ^il nas been added to 
the Hamiltonian H a for a = L [sec Eq. ((2J)] at the start 
of the TD simulation. We used longer leads (-ZVl/r=220) 
and limited the time of the simulations to 100 fs, which 
is sufficient to observe the time propagation before the 
reflections from the finite boundaries start to affect the 
dynamics. The time-dependence of the charge density, 
current and on-site potential, obtained from the closed- 
boundary simulation, is presented in Fig. 01 In contrast 
to our open-boundary simulations, we indeed observed 
qualitatively the same transient regime as in Ref. [13 ■ 
This is mainly characterized by an earlier onset of the 
CB oscillations for larger source-drain voltages and by 
the increase of the width of the first pulse in the time- 
dependence of the Kohn-Sham potential with increasing 
V sd . 



2. Steady-state transport 

In the previous section we demonstrated that, within 
a certain range of parameters, the derivative disconti- 
nuity prevents the quantum dot to evolve towards the 
steady-state. Outside this range, however, a steady-state 
is achievable. Here we determine the steady-state cur- 
rent through the dot for various gate voltages and map 
out the corresponding I- V curves. For situations, where 
the dot is trapped in oscillations, we take as steady-state 
current its time-average in the long-time limit. 

The linear response conductance as a function of V g 
is depicted in Fig. \E\ This is calculated as the finite- 
difference ratio A/o/AT4d close to zero bias (for a very 
low but finite bias AV^d = 0.01 eV) and represents an 
approximation to the zero-bias differential conductance. 
In the non-interacting case, the conductance is composed 
of a single peak centered around V g =1.5 eV, which cor- 
responds to the Fermi level of the leads. This is expected 
from the steady-state picture of transport through a non- 
interacting resonant level. In principle the width of the 
resonance peak is given by the dot-lead hopping integral 
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FIG. 4: (Color online) Real-time evolution of the quantum 
dot in the closed-boundary setup: (a) Charge density of the 
dot (no) for three different values of the source-drain voltage, 
V s d=1.2, 1.3, and 1.4 eV. Current through the dot, 7o, [thick 
lines] and the corresponding Kohn-Sham potential, vks, [thin 
lines] for (b) V sd =1.2 eV, (c) Kd=l-3 eV, and (d) V Bd =lA eV. 
The parameters are the same as in Fig. [3] The source-drain 
voltage is applied as a rigid shift of the on-site energies in the 
left lead. 



7 C . In our TD MPB calculations, however, there is an ad- 
ditional resonance broadening factor (ta) related to the 
dephasing condition in the equations of motion. Its cor- 
responding energy unit, A = /i/ta, can be associated to 
a fictitious temperature, smearing the electronic energy 
distributions in the leads ^ As a result, a suppression of 
the transmission resonance proportional to 1/A is also 
expected. This is the reason of why the amplitude of 
non-interacting resonance conductance in Fig. [5] is below 
one quantum of conductance, Go = 2e 2 /h. 

In the interacting case the Anderson impurity model 
predicts two distinct Coulomb peaks*^ in the conduc- 
tance as a function of the gate voltage^. These are 
manifestation of charge quantization at the dot and corre- 
spond to each of the two integer electron number states, 
in which the dot is inhabited by one or two electrons, 
respectively. Although the ABALDA potential succeeds 
in describing some important properties of strongly cor- 
related systems^ due to the presence of the derivative 
discontinuity, it is a single-particle potential and, as such, 
cannot describe fully these charge states. As a result, the 
gate-voltage dependence of the conductance, calculated 
using the full discontinuous effective potential (wks); does 
not show two distinct peaks. However, it presents a 
structure, bearing the signature of two broadened and 
overlapping peaks (see Fig. [5|. The distance between 
these quasi-peaks increases with increasing U and cor- 
responds to the value of the jump of the on-site poten- 
tial vks[?"Jo] at the derivative discontinuity. In the case 
of the Hartree approximation, the two-peak structure is 
less pronounced and the two resonances merge into an 




V (eV) 

FIG. 5: (Color online) Differential conductance of the dot as 
a function of the gate voltage (V s ) for the Kohn-Sham po- 
tential, dks, [thick solid lines] and for the Hartree potential, 
vh, [dashed lines] with U=l, 2, and 3 eV, and for the non- 
interacting case (thin solid line). The inset shows a compar- 
ison between the density-dependence of vks (solid lines) and 
vh (dashed lines) for the same values of U. Parameters are 
the same as those of Fig. [3] and V s d = 0.01 eV. 



asymmetric plateau. The width of this plateau is also 
proportional to U. 

It should be mentioned that for the TD calculations 
with «ks and for values of V g between the position of the 
U = resonance level V" ros = Ef and V^s ~ U (roughly 
corresponding to the region between the two quasi-pcaks) 
no steady-state is achieved. Hence, the conductance 
curves in this region of V g carry some degree of arbitrari- 
ness, associated with the interpretation of the average TD 
current. In fact, for those gate voltages driving a charge 
density at the dot close to unity, even the calculation of 
the ground-state is problematic from a numerical view- 
point, because of the derivative discontinuity. In such 
cases we used the following iterative procedure. Let V g 
be the value of the gate voltage, for which the ground- 
state (initial) density is calculated self-consistently, while 
V g +5V g is the value of the gate voltage for which the self- 
consistent calculation does not converge. In this case, the 
final density, obtained at the end of the time-dependent 
simulation with V g =V® , is taken as initial density for the 
simulation with V g =V® + SV g . 

In the same way, from the time-averages in the long 
time-limit, we map out the I-V characteristics of the in- 
teracting dot (vks) at a given V g (see Fig. [6]). A remark- 
able feature of the I- V curves is the drop of the current 
(NDC) at large source-drain voltages, which is almost 
negligible for small U but increases with increasing U. 

For all values of U the current initially increases with 
increasing V s d as the dot is charging. It then reaches its 
maximum value as the charge density approaches no = 
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FIG. 6: (Color online) Current through the dot, Io, as a func- 
tion of the source-drain voltage, V s d, for vks and different 
values of U. The horizontal dashed lines represent the cor- 
responding saturation currents Is (see text for the exact def- 
inition). The following parameters are used: 70=— 3.88 eV, 
7c=-0.5 eV, e F =1.5 eV, V g =2.0 eV. 




FIG. 7: Current (a), density (b) and on-site potential (c) of 
the dot as a function of the source-drain voltage, V B d, for vks 
with (7 = 5 eV and for vn with [7 = 9 eV. The inset shows wks 
and vh as functions of the dot density for the corresponding 
values of U. The parameters are the same as those in Fig. [6] 



1. This point corresponds to a threshold source-drain 
voltage V^J, which is roughly the same for all values of U. 
Beyond V^, the system is driven into a dynamical state 
(where the steady-state current is calculated by averaging 
out the oscillations). In the limit of very large V s d, the 
dot recovers its long-time tendency to a steady state and 



the average current saturates. At saturation and beyond 
the dot occupation is above 1 and the on-site potential 
assumes a value above the discontinuity. Hence, the on- 
site energy at the dot is proportional to the the jump of 
the vks at no = 1, i.e. it is proportional to U. 



As discussed in Section MI Al for the non-interacting 
case, the saturation current decreases with increasing the 
dot on-site potential, because of the finite bandwidth of 
the electrodes. For the same reason here the drop of the 
current becomes larger when U increases. In fact a large 
U corresponds to a large value of the steady-state on-site 
potential, which then approaches the electrodes' band- 
edge. In order to confirm this conjecture, we compare the 
saturation current 1$ calculated at finite U, with that for 
[7=0 and V g equal to the steady-state on-site potential 
corresponding to that obtained at the same U. Indeed Is 
matches quite well the value of the current obtained at 
large source-drain voltages in the I- V characteristics of 
the interacting dot (see horizontal dashed lines next to 
each curve in Fig. EJ) - This argument can obviously be re- 
versed, i.e. the NDC cannot be observed, if the variation 
of the on-site potential at the derivative discontinuity, 
determined by U, is much smaller than the electrodes' 
bandwidth. For instance, for the same set of parameters 
used before for the dot+electrodes system, such NDC- 
free situation is found for U = 2 eV (U <C 4|7o| for 
70 = 3.88 eV). In this case the drop of the current above 
is practically negligible. 



Importantly, the NDC displayed in Fig. \6\ is not found 
in I-V's calculated within the Hartrec approximation, 
even for large values of U (see Fig. [7]). When comparing 
calculations at the Hartree level with those performed 
with the complete Kohn-Sham potential we intention- 
ally use different U. These are selected in such a way 
that the value of the potential at hq = 1 is identical in 
the two calculations [see the inset in Fig. Etc)], i.e. in 
such a way that the two calculations give the same sat- 
uration current. At variance with the complete Kohn- 
Sham case, in the Hartree only problem the current, as 
well as the density and the on-site potential, monoton- 
ically increase with V^d until the saturation is reached. 
Based on these numerical results we can argue that the 
sclf-interaction-free shape of the on-site potential ^ks a t 
the dot is a necessary condition for the occurrence of 
the NDC in the I- V. The shallow increase of the on-site 
potential with the charging, produced by the opening 
of the bias window, keeps the resonant level away from 
the electrodes band edge and allows the current to rise. 
Once the on-site charge exceeds no = 1 and the reso- 
nant level energy shoots up towards the band-edge, the 
currents drops. The averaged dynamical current mono- 
tonically approaches its saturation value corresponding 
to a steady-state solution. 
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IV. CONCLUSIONS 

We have investigated the electronic transport through 
a strongly-correlated quantum dot by using a re- 
cently proposed multiple-probe battery method for time- 
dependent simulations of open systems. Our aim was 
two-fold. Firstly, wc wanted to assess the outcomes of 
a TD transport scheme conceptually different from what 
used so far in literature, for a problem involving strong 
electron correlation as in Coulomb blockade. Clearly our 
MPB-based simulations agree well with previous find- 
ings^ In particular we have demonstrated self-sustained 
oscillations in the current, density and effective on-site 
potential, originating from the derivative discontinuity 
of the approximate exchange-correlation potential used. 

As a further aspect we have addressed the question of 
whether the peculiar dynamics obtained from the time- 
dependent simulations can be related to the more acces- 
sible steady-state picture of transport. In particular, we 
have shown the presence of Coulomb peaks in the linear 
response differential conductance and extracted the TD 



version of I- V characteristics, based on the time-averaged 
current through the dot in the long-time limit. The re- 
sulting I-V curves, at a critical voltage, exhibit a drop 
in the average current through the dot. This drop corre- 
sponds to the range of parameters where no steady state 
is found and the dot is in the oscillatory Coulomb block- 
ade state. Such an NDC is however present only when 
the calculation is performed at a DFT level in which the 
potential includes the derivative discontinuity at unitary 
occupation. 
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